Redukcija dimenzionalnosti

In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import shapiro, pearsonr, boxcox
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import r2_score
from sklearn.metrics import confusion_matrix
from scipy.spatial.distance import mahalanobis
from scipy.stats import chi2
from sklearn.preprocessing import MinMaxScaler
from sklearn.neighbors import NearestNeighbors
from sklearn.decomposition import PCA, FactorAnalysis
from sklearn.manifold import TSNE
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
import seaborn as sns

Učitavanje podataka o novčanicama

Izvor: Kaggle

In [2]:
data = pd.read_csv('banknotes.csv')
In [3]:
data.head()
Out[3]:
conterfeit Length Left Right Bottom Top Diagonal
0 0 214.8 131.0 131.1 9.0 9.7 141.0
1 0 214.6 129.7 129.7 8.1 9.5 141.7
2 0 214.8 129.7 129.7 8.7 9.6 142.2
3 0 214.8 129.7 129.6 7.5 10.4 142.0
4 0 215.0 129.6 129.7 10.4 7.7 141.8
In [4]:
data.shape
Out[4]:
(200, 7)
In [5]:
data.describe()
Out[5]:
conterfeit Length Left Right Bottom Top Diagonal
count 200.000000 200.000000 200.000000 200.000000 200.000000 200.000000 200.000000
mean 0.500000 214.896000 130.121500 129.956500 9.417500 10.650500 140.483500
std 0.501255 0.376554 0.361026 0.404072 1.444603 0.802947 1.152266
min 0.000000 213.800000 129.000000 129.000000 7.200000 7.700000 137.800000
25% 0.000000 214.600000 129.900000 129.700000 8.200000 10.100000 139.500000
50% 0.500000 214.900000 130.200000 130.000000 9.100000 10.600000 140.450000
75% 1.000000 215.100000 130.400000 130.225000 10.600000 11.200000 141.500000
max 1.000000 216.300000 131.000000 131.100000 12.700000 12.300000 142.400000

Razdvajanje regularnih novčanica od falsifikata

In [6]:
regular_notes = data[data.conterfeit == 0]
counterfeit_notes = data[data.conterfeit == 1]
In [7]:
print(f'Regular notes: {regular_notes.shape}')
print(f'Counterfeit notes: {counterfeit_notes.shape}')
Regular notes: (100, 7)
Counterfeit notes: (100, 7)

Preprocesiranje

Analiza raspodele atributa, provera normalnosti i, eventualna, transformacija do normalne raspodele kod regularnih novčanica

Atribut: Length

Dužina novčanice

In [8]:
_ = plt.hist(regular_notes.Length)
In [9]:
W, p = shapiro(regular_notes.Length)

# 0.05 -> Kriticna vrednost, za vrednosti manje od 0.05 raspodela nije normalna, sa verovatnocom 0.05.
print(f'Normal distribution? {p > 0.05}, (p={p})')
Normal distribution? True, (p=0.25671127438545227)

Atribute: Left

Duzina leve strane novčanice

In [10]:
_ = plt.hist(regular_notes.Left)
In [11]:
W, p = shapiro(regular_notes.Left)

# 0.05 -> Kriticna vrednost, za vrednosti manje od 0.05 raspodela nije normalna, sa verovatnocom 0.05.
print(f'Normal distribution? {p > 0.05}, (p={p})')
Normal distribution? False, (p=0.00825295876711607)
In [12]:
# Logaritamska transforacija
_, p_log = shapiro(np.log(regular_notes.Left + 1))
print(f'[x -> log(x + 1)]: Normal distribution? {p_log > 0.05}, (p={p_log})\n')

# Transformacija korenovanjem (stepenovanjem)
_, p_sqrt = shapiro(np.sqrt(regular_notes.Left + 1))
print(f'[x -> sqrt(x)]: Normal distribution? {p_sqrt > 0.05}, (p={p_sqrt})\n')

# BOXCOX transformacija
Left_boxcox, Left_lambda = boxcox(regular_notes.Left)
_, p_boxcox = shapiro(Left_boxcox)
print(f'[x -> boxcox(x)]: Normal distribution? {p_boxcox > 0.05}, (p={p_boxcox}, lambda={Left_lambda})')
print(f'Vrednost lambda = {Left_lambda} je prevelika')
[x -> log(x + 1)]: Normal distribution? False, (p=0.008864431641995907)

[x -> sqrt(x)]: Normal distribution? False, (p=0.008556782267987728)

[x -> boxcox(x)]: Normal distribution? True, (p=1.0, lambda=-303.8736987713302)
Vrednost lambda = -303.8736987713302 je prevelika
/usr/local/lib/python3.8/site-packages/scipy/stats/morestats.py:908: RuntimeWarning: divide by zero encountered in log
  return (lmb - 1) * np.sum(logdata, axis=0) - N/2 * np.log(variance)
/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2522: RuntimeWarning: invalid value encountered in double_scalars
  w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom
/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2115: RuntimeWarning: invalid value encountered in double_scalars
  tmp1 = (x - w) * (fx - fv)
/usr/local/lib/python3.8/site-packages/scipy/optimize/optimize.py:2116: RuntimeWarning: invalid value encountered in double_scalars
  tmp2 = (x - v) * (fx - fw)
/usr/local/lib/python3.8/site-packages/scipy/stats/morestats.py:1678: UserWarning: Input data for shapiro has range zero. The results may not be accurate.
  warnings.warn("Input data for shapiro has range zero. The results "

Atribut: Right

Duzina desne strane novčanice

In [13]:
_ = plt.hist(regular_notes.Right)
In [14]:
W, p = shapiro(regular_notes.Right)

# 0.05 -> Kriticna vrednost, za vrednosti manje od 0.05 raspodela nije normalna, sa verovatnocom 0.05.
print(f'Normal distribution? {p > 0.05}, (p={p})')
Normal distribution? False, (p=0.01173424068838358)
In [15]:
# Logaritamska transforacija
_, p_log = shapiro(np.log(regular_notes.Right + 1))
print(f'[x -> log(x + 1)]: Normal distribution? {p_log > 0.05}, (p={p_log})\n')

# Transformacija korenovanjem (stepenovanjem)
_, p_sqrt = shapiro(np.sqrt(regular_notes.Right + 1))
print(f'[x -> sqrt(x)]: Normal distribution? {p_sqrt > 0.05}, (p={p_sqrt})\n')

# BOXCOX transformacija
Right_boxcox, Right_lambda = boxcox(regular_notes.Right)
_, p_boxcox = shapiro(Right_boxcox)
print(f'[x -> boxcox(x)]: Normal distribution? {p_boxcox > 0.05}, (p={p_boxcox}, lambda={Right_lambda})')
print(f'Vrednost lambda = {Right_lambda} je prevelika')
[x -> log(x + 1)]: Normal distribution? False, (p=0.013035133481025696)

[x -> sqrt(x)]: Normal distribution? False, (p=0.012379295192658901)

[x -> boxcox(x)]: Normal distribution? True, (p=1.0, lambda=-291.23915890047374)
Vrednost lambda = -291.23915890047374 je prevelika

Atribut: Bottom

Širina donje margine

In [16]:
_ = plt.hist(regular_notes.Bottom)
In [17]:
W, p = shapiro(regular_notes.Bottom)

# 0.05 -> Kriticna vrednost, za vrednosti manje od 0.05 raspodela nije normalna, sa verovatnocom 0.05.
print(f'Normal distribution? {p > 0.05}, (p={p})')
Normal distribution? False, (p=0.040356215089559555)
In [18]:
# Logaritamska transforacija
_, p_log = shapiro(np.log(regular_notes.Bottom + 1))
print(f'[x -> log(x + 1)]: Normal distribution? {p_log > 0.05}, (p={p_log}) <-\n')

# Transformacija korenovanjem (stepenovanjem)
_, p_sqrt = shapiro(np.sqrt(regular_notes.Bottom + 1))
print(f'[x -> sqrt(x)]: Normal distribution? {p_sqrt > 0.05}, (p={p_sqrt})\n')

# BOXCOX transformacija
Bottom_boxcox, Bottom_lambda = boxcox(regular_notes.Bottom)
_, p_boxcox = shapiro(Bottom_boxcox)
print(f'[x -> boxcox(x)]: Normal distribution? {p_boxcox > 0.05}, (p={p_boxcox}, lambda={Bottom_lambda})')
[x -> log(x + 1)]: Normal distribution? True, (p=0.11221223324537277) <-

[x -> sqrt(x)]: Normal distribution? True, (p=0.07275495678186417)

[x -> boxcox(x)]: Normal distribution? True, (p=0.17343851923942566, lambda=-1.2196806024003763)

Atribut: Top

Širina gornje margine

In [19]:
_ = plt.hist(regular_notes.Top)
In [20]:
W, p = shapiro(regular_notes.Top)

# 0.05 -> Kriticna vrednost, za vrednosti manje od 0.05 raspodela nije normalna, sa verovatnocom 0.05.
print(f'Normal distribution? {p > 0.05}, (p={p})')
Normal distribution? False, (p=0.04984739050269127)
In [21]:
# Logaritamska transforacija
_, p_log = shapiro(np.log(regular_notes.Top + 1))
print(f'[x -> log(x + 1)]: Normal distribution? {p_log > 0.05}, (p={p_log})\n')

# Transformacija korenovanjem (stepenovanjem)
_, p_sqrt = shapiro(np.sqrt(regular_notes.Top + 1))
print(f'[x -> sqrt(x)]: Normal distribution? {p_sqrt > 0.05}, (p={p_sqrt})\n')

# BOXCOX transformacija
Top_boxcox, Top_lambda = boxcox(regular_notes.Top)
_, p_boxcox = shapiro(Top_boxcox)
print(f'[x -> boxcox(x)]: Normal distribution? {p_boxcox > 0.05}, (p={p_boxcox}, lambda={Top_lambda}) <--')
[x -> log(x + 1)]: Normal distribution? False, (p=0.002854172606021166)

[x -> sqrt(x)]: Normal distribution? False, (p=0.013121336698532104)

[x -> boxcox(x)]: Normal distribution? True, (p=0.8009419441223145, lambda=3.4122960497938575) <--

Atribut: Diagonal

Dužina diagonale novčanice

In [22]:
_ = plt.hist(regular_notes.Diagonal)
In [23]:
W, p = shapiro(regular_notes.Diagonal)

# 0.05 -> Kriticna vrednost, za vrednosti manje od 0.05 raspodela nije normalna, sa verovatnocom 0.05.
print(f'Normal distribution? {p > 0.05}, (p={p})')
Normal distribution? False, (p=0.0032950430177152157)
In [24]:
# Logaritamska transforacija
_, p_log = shapiro(np.log(regular_notes.Diagonal + 1))
print(f'[x -> log(x + 1)]: Normal distribution? {p_log > 0.05}, (p={p_log})\n')

# Transformacija korenovanjem (stepenovanjem)
_, p_sqrt = shapiro(np.sqrt(regular_notes.Diagonal + 1))
print(f'[x -> sqrt(x)]: Normal distribution? {p_sqrt > 0.05}, (p={p_sqrt})\n')

# BOXCOX transformacija
Diagonal_boxcox, Diagonal_lambda = boxcox(regular_notes.Diagonal)
_, p_boxcox = shapiro(Diagonal_boxcox)
print(f'[x -> boxcox(x)]: Normal distribution? {p_boxcox > 0.05}, (p={p_boxcox}, lambda={Diagonal_lambda}) <--')
[x -> log(x + 1)]: Normal distribution? False, (p=0.0028021091129630804)

[x -> sqrt(x)]: Normal distribution? False, (p=0.0030389565508812666)

[x -> boxcox(x)]: Normal distribution? True, (p=1.0, lambda=64.37147740204965) <--
/usr/local/lib/python3.8/site-packages/numpy/core/_methods.py:195: RuntimeWarning: overflow encountered in multiply
  x = um.multiply(x, x, out=x)

Pregled raspodela atributa kod regularnih novcanica i falsifikata

Length:

In [25]:
_ = plt.hist(regular_notes.Length)
_ = plt.hist(counterfeit_notes.Length)

Left:

In [26]:
_ = plt.hist(regular_notes.Left)
_ = plt.hist(counterfeit_notes.Left)

Right:

In [27]:
_ = plt.hist(regular_notes.Right)
_ = plt.hist(counterfeit_notes.Right)

log Bottom:

In [28]:
_ = plt.hist(np.log(regular_notes.Right + 1))
_ = plt.hist(np.log(counterfeit_notes.Right + 1))

boxcox Top:

In [29]:
_ = plt.hist(boxcox(regular_notes.Top, lmbda=Top_lambda))
_ = plt.hist(boxcox(counterfeit_notes.Top, lmbda=Top_lambda))
In [30]:
_ = plt.hist(boxcox(regular_notes.Diagonal, lmbda=Diagonal_lambda))
_ = plt.hist(boxcox(counterfeit_notes.Diagonal, lmbda=Diagonal_lambda))

Z Score

Napomena: Promenom vrednosti praga sa 3 na niže vrednosti (2.5) dobija se bolji rezultat predikcije (probati). Takodje, korišćenjem netransformisane vrednosti diagonale dobijaju se bolji rezultati jer je raspodela vrlo bliska normalnoj (probati). Ipak, kako formalni uslov normalnosti nije ispunjen, korišćene su vrednosti dobijene BOXCOX transformacijom.

In [31]:
bc_Diagonal = boxcox(regular_notes.Diagonal, lmbda=Diagonal_lambda)
mu = bc_Diagonal.mean()
sigma = bc_Diagonal.std()

print(f'Mean = {mu}, Sigma = {sigma}')

y = [] # Stvarne oznake klase
y_pred = [] # Predvidjene oznake klase

for index, row in data.iterrows():
    # Po definiciji BOXCOX transformacije
    x = (row.Diagonal**Diagonal_lambda - 1) / Diagonal_lambda

    y.append(row.conterfeit)
    
    z = abs(x - mu) / sigma
    
    if z > 3:
        y_pred.append(1)
    else:
        y_pred.append(0)
Mean = 4.4707410261888524e+136, Sigma = 8.550314200577817e+135
In [32]:
confusion_matrix(y, y_pred)
Out[32]:
array([[99,  1],
       [39, 61]])
In [33]:
print(f'Score: {r2_score(y, y_pred)}')
Score: 0.19999999999999996

Mahalanobisovo rastojanje

In [34]:
X = regular_notes.loc[:,['Length','Diagonal']]
X.Diagonal = boxcox(X.Diagonal, lmbda=Diagonal_lambda)

y = []
y_pred = []

mu = X.mean()
    
V = X.cov()
VI = np.linalg.pinv(V)

p = 0.95
df = X.shape[1]

critical_value = chi2.ppf(p, df)
print(f'Critical value: {critical_value}')

for index, row in data.iterrows():
    x = [
        row.Left, 
        (row.Diagonal**Diagonal_lambda - 1) / Diagonal_lambda]
    y.append(row.conterfeit)
    
    M = mahalanobis(x, mu, VI)
    
    if M**2 > critical_value:
        y_pred.append(1)
    else:
        y_pred.append(0)
Critical value: 5.991464547107979
In [35]:
confusion_matrix(y, y_pred)
Out[35]:
array([[99,  1],
       [ 8, 92]])
In [36]:
print(f'Score: {r2_score(y, y_pred)}')
Score: 0.8200000000000001

t-SNE

t-SNE se može koristiti samo za vizualizaciju

In [37]:
# Izdvajanje atributa koji nisu oznaka klase
data_attributes = data.iloc[:, 1:] # nulta kolona je oznaka klase
data_attributes.head()
Out[37]:
Length Left Right Bottom Top Diagonal
0 214.8 131.0 131.1 9.0 9.7 141.0
1 214.6 129.7 129.7 8.1 9.5 141.7
2 214.8 129.7 129.7 8.7 9.6 142.2
3 214.8 129.7 129.6 7.5 10.4 142.0
4 215.0 129.6 129.7 10.4 7.7 141.8
In [38]:
tsne = TSNE(n_components=2, perplexity=30, n_iter=20000)
tsne_data_transformed = tsne.fit_transform(data_attributes)
In [39]:
X = tsne_data_transformed[:,0]
Y = tsne_data_transformed[:,1]

_ = plt.scatter(X, Y, c=data.conterfeit)

PCA - Principal Component Analysis

In [40]:
pca = PCA()
pca_data_transformed = pca.fit_transform(data_attributes)

Analiza udela objasnjene varijanse pomoću svake od kompomenti:

In [41]:
evr = pca.explained_variance_ratio_

for i in range(evr.shape[0]):
    print(f'PC{i + 1}: {evr[i] * 100}%')

_ = plt.bar([i + 1 for i in range(evr.shape[0])], evr)
PC1: 66.75169871502327%
PC2: 20.81597102092934%
PC3: 5.415542450842821%
PC4: 4.330827173088247%
PC5: 1.8958192790454902%
PC6: 0.7901413610708573%

Pomoću prve dve komponente može se objasniti ~87% udela varijanse

In [42]:
X = pca_data_transformed[:,0]
Y = pca_data_transformed[:,1]

_ = plt.scatter(X, Y, c=data.conterfeit)

Izračunavanje komponenti vršiće se na osnovu regularnih novčanica, kako falsifikati ne bi imali uticaj u proceni parametara raspodele. Ovo će biti vršeno na navedeni način samo u ovom specifičnom slučaju prepoznavanja elemenata van granica, u slučaju klasifikacije biće transformisan ceo skup podataka.

In [43]:
pca_regular_data_transformed = pca.fit_transform(regular_notes.iloc[:,1:])

Pomocu prve tri komponente može se objasniti 85% udela varijanse

In [44]:
evr = pca.explained_variance_ratio_

for i in range(evr.shape[0]):
    print(f'PC{i + 1}: {evr[i] * 100}%')

_ = plt.bar([i + 1 for i in range(evr.shape[0])], evr)
PC1: 47.741696602756825%
PC2: 24.8928965668078%
PC3: 12.860106340416632%
PC4: 6.044208478138688%
PC5: 5.555616868113948%
PC6: 2.905475143766102%

Provera normalnosti prve tri komponente:

In [45]:
_ = plt.hist(pca_regular_data_transformed[:,0])
_, p = shapiro(pca_regular_data_transformed[:,0])
print(f'Normal distribution? {p > 0.05}')
Normal distribution? False
In [46]:
log_const = 5 # Konstanta, kako vrednosti ne bi bile negativne za potrebe logaritma
plt.hist(np.log(pca_regular_data_transformed[:,0] + log_const))
_, p_log = shapiro(np.log(pca_regular_data_transformed[:,0] + log_const))
print(f'PC1: Normal distribution? {p_log > 0.05}')
PC1: Normal distribution? True
In [47]:
_ = plt.hist(pca_regular_data_transformed[:,1])
_, p = shapiro(pca_regular_data_transformed[:,1])
print(f'PC2: Normal distribution? {p > 0.05}')
PC2: Normal distribution? False
In [48]:
log_const = 5 # Konstanta, kako vrednosti ne bi bile negativne za potrebe logaritma
# Napomena: Konstanta ima uticaj na oblik raspodele tako da treba biti pazljiv sa vrednostima
# sto manja vrednost - to manji uticaj

plt.hist(np.log(pca_regular_data_transformed[:,1] + log_const))
_, p_log = shapiro(np.log(pca_regular_data_transformed[:,1] + log_const))
print(f'Normal distribution? {p_log > 0.05}')
Normal distribution? True
In [49]:
_ = plt.hist(pca_regular_data_transformed[:,2])
_, p = shapiro(pca_regular_data_transformed[:,2])
print(f'PC2: Normal distribution? {p > 0.05}')
PC2: Normal distribution? True
In [50]:
pca_all_data_transformed = pca.transform(data_attributes)
pca_all_data_transformed_df = pd.DataFrame(data=pca_all_data_transformed[:,:3], columns=['PC1','PC2','PC3'])

# Transformacije to normalne raspodele
pca_all_data_transformed[:,0] = np.log(pca_all_data_transformed[:,0] + log_const)
pca_all_data_transformed[:,1] = np.log(pca_all_data_transformed[:,1] + log_const)
In [51]:
X = pd.DataFrame(data=pca_all_data_transformed[:100,:3], columns=['PC1','PC2','PC3'])
X.head()
Out[51]:
PC1 PC2 PC3
0 1.760318 1.850750 0.247163
1 1.671717 1.468901 0.075238
2 1.746456 1.501675 -0.383679
3 1.460607 1.476851 -0.383502
4 2.107545 1.520611 0.192872
In [52]:
y = []
y_pred = []

mu = X.mean()
    
V = X.cov()
VI = np.linalg.pinv(V)

p = 0.95
df = X.shape[1]

critical_value = chi2.ppf(p, df)
print(f'Critical value: {critical_value}')

for index, row in pca_all_data_transformed_df.iterrows():
    x = row
    y.append(int(data.iloc[index,:].conterfeit))
    
    M = mahalanobis(x, mu, VI)
    
    if M**2 > critical_value:
        y_pred.append(1)
    else:
        y_pred.append(0)
Critical value: 7.814727903251179
In [53]:
confusion_matrix(y, y_pred)
Out[53]:
array([[ 96,   4],
       [  0, 100]])
In [54]:
r2_score(y, y_pred)
Out[54]:
0.92

Redukcija za potrebe klasifikacije

LDA - Linear Discriminant Analysis

In [55]:
lda = LinearDiscriminantAnalysis()
lda_data_transformed = lda.fit_transform(data_attributes, data.conterfeit) # atributi, oznake klase
In [56]:
lda_data_transformed
_ = plt.scatter(lda_data_transformed, [0 for i in range(lda_data_transformed.shape[0])], c=data.conterfeit) # Podaci su jednodimenzioni

U ovom, trivijalnom slučaju, klasifikacija se može vršiti jednostavnim stablom dubine 1 (jedno pitanje):

In [57]:
y = []
y_pred = []

for i in range(lda_data_transformed.shape[0]):
    x = lda_data_transformed[i]
    
    if i < 100:
        y.append(0)
    else:
        y.append(1)
    
    if x[0] < 0:
        y_pred.append(0)
    else:
        y_pred.append(1)
In [58]:
confusion_matrix(y, y_pred)
Out[58]:
array([[ 99,   1],
       [  0, 100]])
In [59]:
r2_score(y, y_pred)
Out[59]:
0.98
In [60]:
num_components = data_attributes.shape[1]
mdl = FactorAnalysis(n_components=num_components)
fa_data_transformed = mdl.fit_transform(data_attributes)
In [61]:
fa_data_transformed = pd.DataFrame(data = fa_data_transformed, columns=[f'Factor {i + 1}' for i in range(num_components)])
fa_data_transformed.head()
Out[61]:
Factor 1 Factor 2 Factor 3 Factor 4 Factor 5 Factor 6
0 -0.544707 1.105304 0.778478 0.0 0.0 0.0
1 1.118511 -0.390437 0.114523 0.0 0.0 0.0
2 1.114295 -0.332905 0.387341 0.0 0.0 0.0
3 1.271035 0.020936 -0.229645 0.0 0.0 0.0
4 0.929099 -0.986241 1.420074 0.0 0.0 0.0

Može se uočiti da je dovoljno sačuvati samo prva tri faktora

Koeficijenti faktora u atributima:

In [62]:
pd.DataFrame(data=mdl.components_, columns=data_attributes.columns).head()
Out[62]:
Length Left Right Bottom Top Diagonal
0 0.033030 -0.231974 -0.275890 -1.040123 -0.394726 0.907176
1 0.171514 0.112144 0.103826 -0.366359 0.135212 0.081033
2 0.066429 0.031747 0.034921 0.257121 -0.264109 0.174766
3 -0.000000 0.000000 0.000000 -0.000000 -0.000000 0.000000
4 -0.000000 0.000000 -0.000000 -0.000000 -0.000000 -0.000000

Uporediti dobijene koeficijente sa korelacijama medju atributima:

In [63]:
_ = sns.heatmap(data_attributes.corr())